Maps
## Saving 16 x 9.89 in image
## Simple feature collection with 6 features and 19 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -171.0125 ymin: 57.8875 xmax: -101.6625 ymax: 65.5875
## Geodetic CRS: WGS 84
## # A tibble: 6 × 20
## City coordinates `Unnamed: 0` mean_ndvi mean_surface_water mean_uhi_DayNight
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Broch… [-101.6624… 202 1866. 3 NA
## 2 Eagle [-141.1958… 305 2408. 2.5 NA
## 3 Emmon… [-164.5208… 373 2971. NA NA
## 4 Lavre… [-171.0124… 504 794. 2.33 NA
## 5 Lutse… [-110.7375… 91 1330. 2.75 NA
## 6 Mekor… [-166.1958… 337 5845 2 NA
## # ℹ 14 more variables: mean_uhi_daytime <dbl>, mean_uhi_Nighttime <dbl>,
## # Continent <chr>, Class <chr>, Skewness <dbl>, Dip_Statistic <dbl>,
## # Mean_Height <dbl>, shape <fct>, shape_emoji <chr>, shape_color <chr>,
## # size_emoji <dbl>, lon <dbl>, lat <dbl>, geometry <POINT [°]>
## Warning: There were 2 warnings in `stopifnot()`.
## The first warning was:
## ℹ In argument: `logMean = case_when(...)`.
## Caused by warning:
## ! NaNs produced
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## `summarise()` has grouped output by 'var1'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'var1'. You can override using the
## `.groups` argument.
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -158.0104 ymin: -45.88905 xmax: 178.023 ymax: 69.41286
## Geodetic CRS: WGS 84
## # A tibble: 6 × 4
## # Groups: var1 [2]
## var1 Shape Meanvar geometry
## <chr> <chr> <dbl> <MULTIPOINT [°]>
## 1 " Daytime UHI (°C)" Diamond 0.0254 ((-11.73729 7.961042), (15.37292 10.2…
## 2 " Daytime UHI (°C)" Hourglass 0.818 ((-41.7813 -22.36964), (-42.17087 -22…
## 3 " Daytime UHI (°C)" Pyramid 0.803 ((-42.52779 -22.27071), (-43.37084 -2…
## 4 " Nighttime UHI (°C)" Diamond 0.522 ((-11.73729 7.961042), (15.37292 10.2…
## 5 " Nighttime UHI (°C)" Hourglass 0.656 ((-41.7813 -22.36964), (-42.17087 -22…
## 6 " Nighttime UHI (°C)" Pyramid 0.487 ((-42.52779 -22.27071), (-43.37084 -2…
## Saving 16 x 9.89 in image
#Violin plots of variables UHI, NDVI, Surface water
ggplot(data = meanval_log, aes(x=Shape, y=logMean , group=interaction(Shape,var1), fill=var1)) +
geom_violin()+
theme_minimal()+
labs(title = "Violin plot of UHI, NDVI, Surface water mean values " ,
subtitle = "Based on city shapes")+
xlab("City Shape")+
ylab("Log Mean ")+
facet_wrap(. ~ var1 , scales = "free", nrow = 2)+
theme(legend.position = "none",
panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_rect(fill = "azure"),
text = element_text(size = 15, face = "bold") ,
axis.text.x=element_text(angle=30),
plot.title = element_text(size = 25, face = "bold"),
plot.subtitle = element_text(size = 15, face = "bold")
) # _3_4_5
ggsave("../data/11_Output_R/ViolinPlotsOfVariablesLog.png")
## Saving 16 x 9.89 in image
library("ggpubr")
# order = c("ctrl", "trt1", "trt2"),
m1_sum <- meanval_sum %>%
filter( var1=='Both Day and Night UHI (°C)')
m1 <- meanval %>%
filter( var1=='Both Day and Night UHI (°C)')
msel <- dplyr::select(m1, Shape, Mean)
ggline(msel, x = "Shape", y = "Mean",
order = c("Diamond", "Hourglass", "Pyramid"),
add= "mean_se", group = 1,
ylab = "Mean", xlab = "City Shape")+
geom_line(data = m1_sum , aes(x=Shape, y=Meanvar, group=1))+
labs(title = "Anova plot of UHI Day and Night value " ,
subtitle = "Based on city shapes")
ggsave("../data/11_Output_R/Anova.png")
## Saving 16 x 9.89 in image
res.aov <- aov(Mean ~ Shape, data = msel)
# Summary of the analysis
summary(res.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Shape 2 42 21.157 15.82 1.44e-07 ***
## Residuals 3745 5008 1.337
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kruskal.test(Mean ~ Shape, data = msel)
##
## Kruskal-Wallis rank sum test
##
## data: Mean by Shape
## Kruskal-Wallis chi-squared = 35.428, df = 2, p-value = 2.027e-08
pairwise.wilcox.test(msel$Mean, msel$Shape,
p.adjust.method = "BH")
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: msel$Mean and msel$Shape
##
## Diamond Hourglass
## Hourglass 0.0047 -
## Pyramid 0.0060 5e-07
##
## P value adjustment method: BH
m2_sum <- meanval_sum %>%
filter( var1==' Daytime UHI (°C)')
m2 <- meanval %>%
filter( var1==' Daytime UHI (°C)')
msel2 <- dplyr::select(m2, Shape, Mean)
ggline(msel2, x = "Shape", y = "Mean",
order = c("Diamond", "Hourglass", "Pyramid"),
add= "mean_se", group = 1,
ylab = "Mean", xlab = "City Shape")+
geom_line(data = m2_sum , aes(x=Shape, y=Meanvar, group=1))+
labs(title = "Anova plot of UHI Day value " ,
subtitle = "Based on city shapes")
ggsave("../data/11_Output_R/AnovaDay.png")
## Saving 16 x 9.89 in image
res.aov <- aov(Mean ~ Shape, data = msel2)
# Summary of the analysis
summary(res.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Shape 2 9 4.714 4.082 0.0169 *
## Residuals 3745 4325 1.155
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kruskal.test(Mean ~ Shape, data = msel2)
##
## Kruskal-Wallis rank sum test
##
## data: Mean by Shape
## Kruskal-Wallis chi-squared = 9.4101, df = 2, p-value = 0.00905
pairwise.wilcox.test(msel$Mean, msel$Shape,
p.adjust.method = "BH")
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: msel$Mean and msel$Shape
##
## Diamond Hourglass
## Hourglass 0.0047 -
## Pyramid 0.0060 5e-07
##
## P value adjustment method: BH
m3_sum <- meanval_sum %>%
filter( var1==' Nighttime UHI (°C)')
m3 <- meanval %>%
filter( var1==' Nighttime UHI (°C)')
msel3 <- dplyr::select(m3, Shape, Mean)
ggline(msel3, x = "Shape", y = "Mean",
order = c("Diamond", "Hourglass", "Pyramid"),
add= "mean_se", group = 1,
ylab = "Mean", xlab = "City Shape")+
geom_line(data = m3_sum , aes(x=Shape, y=Meanvar, group=1))+
labs(title = "Anova plot of UHI Day value " ,
subtitle = "Based on city shapes")
ggsave("../data/11_Output_R/AnovaNight.png")
## Saving 16 x 9.89 in image
res.aov <- aov(Mean ~ Shape, data = msel2)
# Summary of the analysis
summary(res.aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Shape 2 9 4.714 4.082 0.0169 *
## Residuals 3745 4325 1.155
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kruskal.test(Mean ~ Shape, data = msel2)
##
## Kruskal-Wallis rank sum test
##
## data: Mean by Shape
## Kruskal-Wallis chi-squared = 9.4101, df = 2, p-value = 0.00905
pairwise.wilcox.test(msel$Mean, msel$Shape,
p.adjust.method = "BH")
##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: msel$Mean and msel$Shape
##
## Diamond Hourglass
## Hourglass 0.0047 -
## Pyramid 0.0060 5e-07
##
## P value adjustment method: BH
library(GGally)
library(dplyr)
corval <- cty_sf |>
mutate(
Shape = case_when( Class == "Diamond" ~ 1,
Class == "Hourglass" ~ 2 ,
Class == "Pyramid" ~ 3,
Class == "Inverse Pyramid" ~ 4,
TRUE ~ 0 )
)
corval <- dplyr::select(corval,'mean_uhi_DayNight', 'mean_uhi_daytime',
'mean_uhi_Nighttime', 'mean_ndvi' , 'mean_surface_water',
'Shape' )
corval <- corval %>% drop_na()
ggpairs( corval, columns=c('mean_uhi_DayNight', 'mean_uhi_daytime',
'mean_uhi_Nighttime', 'mean_ndvi' , 'mean_surface_water', 'Shape'))
#ggpairs( cty_sf, columns=c('mean_uhi_DayNight', 'Class', 'mean_ndvi' , 'mean_surface_water')
#)
ggsave("../data/11_Output_R/pairsDayNight.png")
## Saving 16 x 9.89 in image
ggpairs( corval, columns=c('mean_uhi_DayNight', 'Shape'))
#ggpairs( cty_sf, columns=c('mean_uhi_DayNight', 'Class', 'mean_ndvi' , 'mean_surface_water')
#)
ggsave("../data/11_Output_R/pairs.png")
## Saving 16 x 9.89 in image
corval <- dplyr::select(corval,'mean_uhi_DayNight', 'mean_uhi_daytime',
'mean_uhi_Nighttime', 'mean_ndvi' , 'mean_surface_water',
'Shape')
corval <- corval %>% dplyr::select(-geometry)
corval <- corval %>%
mutate( Shp = Shape
)
corval2 <- st_drop_geometry(corval)
#cormsel = subset(cormsel, select = -c(Shape,geometry) )
corval2 <- corval2 %>% dplyr::select(-Shape) # geometry)
cormsel1 <- msel %>% dplyr::select(-geometry)
#cormsel1 = subset(cormsel1, dplyr::select = c(Shape,Mean) )
cormsel1 <- cormsel1 %>%
mutate(
Shp = case_when( Shape == "Diamond" ~ 1,
Shape == "Hourglass" ~ 2 ,
Shape == "Pyramid" ~ 3,
Shape == "Inverse Pyramid" ~ 4,
TRUE ~ 0 )
)
cormsel2 <- st_drop_geometry(cormsel1)
#cormsel = subset(cormsel, select = -c(Shape,geometry) )
cormsel2 <- cormsel2 %>% dplyr::select(-Shape) # geometry)
cormsel2 = subset(cormsel2, select = c(Shp, Mean) )
head(cormsel2)
## # A tibble: 6 × 2
## Shp Mean
## <dbl> <dbl>
## 1 3 -6.50
## 2 3 -5.11
## 3 3 -4.78
## 4 2 -4.25
## 5 2 -4.00
## 6 2 -3.85
#ggpairs( cormsel, columns=c('Mean', 'Shp') )
#ggsave("../data/11_Output_R/pairsDayNight_Shape.png")
library("corrgram")
#head(corval)
#vars2 <- c("mean_uhi_daytime","mean_uhi_Nighttime","mean_ndvi", "mean_surface_water")
vars2 <- c("Shp","Mean")
corrgram( cormsel2, main="UHI and City Shape")
ggsave("../data/11_Output_R/Corrgram_uhi_city.png")
## Saving 16 x 9.89 in image
corrgram( corval2, main="UHI and all other variables")
ggsave("../data/11_Output_R/Corrgram_uhi_all_variables.png")
## Saving 16 x 9.89 in image
## class : SpatRaster
## dimensions : 43, 68, 2 (nrow, ncol, nlyr)
## resolution : 0.002694934, 0.002694934 (x, y)
## extent : 9.849982, 10.03324, 56.98437, 57.10025 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : Aalborg_uhi.tif
## names : Aalborg_uhi_1, Aalborg_uhi_2
## min values : -1.186516, -0.2753552
## max values : 2.113822, 0.9442680
## [1] "Aalborg_uhi_1"
## class : SpatRaster
## dimensions : 52, 82, 1 (nrow, ncol, nlyr)
## resolution : 0.002245788, 0.002245788 (x, y)
## extent : 9.850027, 10.03418, 56.98238, 57.09917 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : Aalborg_NDVI.tif
## name : NDVI_mean
## min value : -1306.000
## max value : 7139.435
## [1] "NDVI"
## class : SpatRaster
## dimensions : 43, 68, 1 (nrow, ncol, nlyr)
## resolution : 0.002694946, 0.002694946 (x, y)
## extent : 9.850027, 10.03328, 56.98463, 57.10051 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : Aalborg_SWATER.tif
## name : waterClass
## min value : 1
## max value : 3
## [1] "SWATER"
## Saving 16 x 9.89 in image
## Saving 16 x 9.89 in image
## Saving 16 x 9.89 in image
## Saving 16 x 9.89 in image
## Saving 16 x 9.89 in image
# Not Used
#Boxplots of variables UHI, NDVI, Surface water
ggplot(data = meanval, aes(x=Shape, y=Mean , group=interaction(Shape,var1), fill=var1)) +
geom_boxplot( width =0 )+
stat_boxplot(geom = "errorbar") +
geom_point(data = meanval_sum , aes(x=Shape, y=Meanvar, size=2))+
geom_line(data = meanval_sum , aes(x=Shape, y=Meanvar, group=var1))+
theme_minimal()+
labs(title = "Boxplot of UHI, NDVI, Surface water global mean values " ,
subtitle = "Based on city shapes")+
xlab("City Shape")+
ylab("Mean ")+
facet_wrap(. ~ var1 , scales = "free", nrow = 2)+
theme(legend.position = "none",
panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_rect(fill = "azure"),
text = element_text(size = 15, face = "bold") ,
axis.text.x=element_text(angle=30),
plot.title = element_text(size = 25, face = "bold"),
plot.subtitle = element_text(size = 15, face = "bold")
) # _3_4_5
ggsave("../data/11_Output_R/BoxPlotsErrorbarOfVariables.png")
## Saving 16 x 9.89 in image
#Boxplots of variables UHI, NDVI, Surface water log values
ggplot(data = meanval_log, aes(x=Shape, y=logMean , group=interaction(Shape,var1), fill=var1)) +
geom_boxplot( width =0 )+
stat_boxplot(geom = "errorbar") +
geom_point(data = meanval_log_sum , aes(x=Shape, y=Meanvar, size=2))+
geom_line(data = meanval_log_sum , aes(x=Shape, y=Meanvar, group=var1))+
theme_minimal()+
labs(title = "Boxplot of UHI, NDVI, Surface water global log mean values " ,
subtitle = "Based on city shapes")+
xlab("City Shape")+
ylab("LogMean ")+
facet_wrap(. ~ var1 , scales = "free", nrow = 2)+
theme(legend.position = "none",
panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_rect(fill = "azure"),
text = element_text(size = 15, face = "bold") ,
axis.text.x=element_text(angle=30),
plot.title = element_text(size = 25, face = "bold"),
plot.subtitle = element_text(size = 15, face = "bold")
) # _3_4_5
ggsave("../data/11_Output_R/BoxPlotsErrorbarOfVariables_log.png")
## Saving 16 x 9.89 in image